USE ElliDef
USE NodeInfoDef
USE metric

USE Elliptic_Solvers
IMPLICIT NONE
TYPE (NodeInfo) :: Info(2)
INTEGER N,error
REAL (KIND=qPrec) :: dx,dy,dz
REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: p,derp,u,v,w,div,pl
REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:,:) :: ut,vt,wt
! first thing allocate structures for sparse matrixes
INTEGER i,j,k,nx,ny,nz
!INTEGER, PARAMETER :: solver=PARDISO_SOLVER
!INTEGER, PARAMETER :: solver=IT_SD
!INTEGER, PARAMETER :: solver=IT_CG
!INTEGER, PARAMETER :: solver=IT_FOM
!INTEGER, PARAMETER :: solver=DCG_SOLVER
INTEGER, PARAMETER :: solver=IT_GMRES
INTEGER, PARAMETER :: how_many_steps=5
LOGICAL iterate







DO i=1,2
! populate Info structures
Info(i)%mX=(/200,40,20/)

nx=Info(i)%mX(1);ny=Info(i)%mX(2);nz=Info(i)%mX(3)

Info(i)%nDim=3 ! 3D problem
! now define geometry of the problem
Info(i)%xLower=(/-15000./2200.,-15000./2200.,0./)

Info(i)%xUpper=(/20000./2200,15000./2200,1./)

Info(i)%dX=(Info(i)%xUpper(1:3)-Info(i)%xLower(1:3))/Info(i)%mX
Info(i)%dmax=0.

! assume all Neumann conditions
Info(i)%mthbc(1:6)=0
END DO

!now copy info to elli structure


ALLOCATE(Info(1)%ElliInfo) ! allocate space for ellikit 
Info(1)%ElliInfo%ndim=Info(1)%ndim ! dimension of problem
Info(1)%ElliInfo%mX=Info(1)%mX ! size of grid
Info(1)%ElliInfo%dX=Info(1)%dX ! grid increments
Info(1)%ElliInfo%mthbc=Info(1)%mthbc ! definition of boundary conditions




   Call setprob(Info(1))
PRINT*,'Dimension of problem',Info(1)%mX

CALL DEFINE_TOPO(Info(1)) ! define the quantities used to calculate the metric

error=CREATE_OPERATORS(Info(1)%ElliInfo)

! Initially zero in the boundary conditions to be safe
DO i=1,2*Info(1)%nDim
   Info(1)%ElliInfo%bc(i)%p=zero

END DO

! check consistency of BCs
IF(.NOT.check(Info(1)%ElliInfo,CHECK_BC)) THEN
   PRINT*,'BC are not consistent'
   STOP
END IF
Info(1)%ElliInfo%qe=zero
! assign local pointers to edge velocities and fill in with random values
u=>Info(1)%ElliInfo%qe(1:nx+1,1:ny,1:nz,1);
v=>Info(1)%ElliInfo%qe(1:nx,1:ny+1,1:nz,2);
w=>Info(1)%ElliInfo%qe(1:nx,1:ny,1:nz+1,3)
ALLOCATE(ut(1:ny,1:nz,2),vt(1:nx,1:nz,2),wt(1:nx,1:ny,2))



CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)
CALL RANDOM_NUMBER(w);
! random values for edge velocities boundary conditions
!!$CALL RANDOM_NUMBER(ut)
!!$CALL RANDOM_NUMBER(vt)
!!$CALL RANDOM_NUMBER(wt)
ut=zero;vt=zero;wt=zero
        

! only need OP= o use IT wth PARDISO precon

i=POISSON_SOLVERS(Info(1)%ElliInfo,INIT,solver,OP=CSR_PRECOND)

Info(1)%ElliInfo%it%ipar_it(2)=0 ! no precond
ALLOCATE(Info(2)%ElliInfo) ! allocate space for ellikit 

Info(2)%ElliInfo%ndim=Info(1)%ndim ! dimension of problem
Info(2)%ElliInfo%mX=Info(1)%mX ! size of grid
Info(2)%ElliInfo%dX=Info(1)%dX ! grid increments
Info(2)%ElliInfo%mthbc=Info(1)%mthbc ! definition of boundary conditions

CALL setprob(Info(2))
CALL DEFINE_TOPO(Info(2)) ! define the quantities used to calculate the metric

error=CREATE_OPERATORS_THIN(Info(2)%ElliInfo)

DEALLOCATE(Info(2)%ELLIINFO%qe);
Info(2)%ElliInfo%qe=>Info(1)%ElliInfo%qe

i=POISSON_SOLVERS_THIN(Info(2)%ElliInfo,INIT)
CALL hodge_thin(Info(2)%ElliInfo,1,.false.,ut,vt,wt)

ITERATE=.true.;i=1
DO WHILE(ITERATE.AND.i<10000)
PRINT*,'iteration ',i
i=i+1


Info(1)%ElliInfo%bc(1)%p(1:ny,1:nz)=u(1,:,:)-ut(:,:,1) ! u*-utarget
Info(1)%ElliInfo%bc(2)%p(1:ny,1:nz)=u(nx+1,:,:)-ut(:,:,2)
Info(1)%ElliInfo%bc(3)%p(1:nx,1:nz)=v(:,1,:)-vt(:,:,1)
Info(1)%ElliInfo%bc(4)%p(1:nx,1:nz)=v(:,ny+1,:)-vt(:,:,2)

Info(1)%ElliInfo%bc(5)%p(1:nx,1:ny)=w(1:nx,1:ny,1)-wt(:,:,1)
Info(1)%ElliInfo%bc(6)%p(1:nx,1:ny)=w(1:nx,1:ny,nz+1)-wt(:,:,2)
!!$


PRINT*,'calculates divergence'
   

CALL take_div(Info(1)%ElliInfo,.true.)


div=>Info(1)%ElliInfo%div
PRINT*,'divergence stat'

!!$


PRINT*,'max divergence before = ',maxval(abs(div))

PRINT*,'rms divergence before = ',sqrt(SUM(div**2)/(nx*ny*nz))
IF(solver==DCG_SOLVER) j=POISSON_SOLVERS(Info(1)%ElliInfo,INIT,solver) ! need to reinit solver
Info(1)%ElliInfo%pold=zero
j=POISSON_SOLVERS(Info(1)%ElliInfo,SOLVE,solver,OP=CSR_PRECOND)

! subtract gradient
CALL project(Info(1)%ElliInfo,.true.)

PRINT*,'calculates divergence'
10 CALL take_div(Info(1)%ElliInfo,.false.)
div=>Info(1)%ElliInfo%div

PRINT*,'max divergence after = ',maxval(abs(div))

PRINT*,'rms divergence after = ',sqrt(SUM(div**2)/(nx*ny*nz))
ITERATE=sqrt(SUM(div**2)/(nx*ny*nz))>1e-10
IF(ITERATE.AND.MOD(i,how_many_steps)==0) CALL hodge_thin(Info(2)%ElliInfo,2,.false.,ut,vt,wt) ! every how_many_steps apply the accelerator
ENDDO

i=POISSON_SOLVERS(Info(1)%ElliInfo,FINISH,solver,OP=CSR_PRECOND)
i=POISSON_SOLVERS_THIN(Info(2)%ElliInfo,FINISH) !wrap up accelerator 
PRINT*,'parison_mopup'
NULLIFY(Info(2)%ElliInfo%qe);ALLOCATE(Info(2)%ElliInfo%qe(1,1,1,1))
DO i=1,2
   error=FREEUP_OPERATORS(Info(i)%ElliInfo)

DEALLOCATE(Info(i)%topo%h,Info(i)%topo%hx,Info(i)%topo%hy,Info(i)%topo%hxx,Info(i)%topo%hyy)
DEALLOCATE(Info(i)%topo)
DEALLOCATE(Info(i)%ElliInfo)
END DO
DEALLOCATE(ut,vt,wt)
STOP

END

